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PREFACE 


This is the final report on the Evaluation of Digital Correction Tech- 
niques for ERTS Images. This report fulfills the deliverable documentation 
requirements for Item 5, Contract NAS5-21814. Additionally, this report 
summarizes work completed since the last progress report. 

During the course of the study, TRW has evaluated various digital 
correction techniques applied to bulk CCT MSS and RBV data in four primary 
areas: 


• processing - induced distortions 

• reseau and ground control point location accuracy 

• geometric distortion correction accuracy 

• throughput 

TRWs CDC 6400 computer was utilized as the major tool in evaluating these 
techniques. The study was not explicitly concerned with image enhancement, 
but rather concentrated on the various elements of ERTS data processing. 

While the various tasks have not been combined to produce a contiguous data 
flow required in an operational facility, the processing modules described 
in this report have been evaluated individually in terms of throughput. 

The principle performer of this study was Dr. S. Rifman, with Mr. J. 
Taber and Dr. K. Simon providing overall direction. Dr. Rifman was period- 
ically supported by M. Benesch, A. Hung, D. McKinnon, J. Rau, W. Spence, 
and K. Ziedman. Ms. Benesch particularly made valuable contributions to the 
RBV effort. Dr. R. Caron provided helpful consultation for the Kalman Filter. 

Deep appreciation is extended toward Mr. Bernard Peavey, the Science 
Monitor, who has provided much valuable support throughout this contract. 
Grateful thanks are due Messrs. Tom Burger and Mel Byerly of the U.S. Geo- 
logical Survey who supplied the necessary precision geodetic ground truth 
data. 
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1 

INTRODUCTION AND 
SUMMARY 

1.1 STUDY SCOPE 

The ultimate value of satellite-collected Earth resources data will 
depend on the speed and accuracy with which meaningful information can be 
extracted. TRW's study objective has been to verify that the prerequisite 
geometric and photometric correction of ERTS-1 MSS and RBV data can be 
speedily and accurately accomplished with all-digital image processing 
techniques. 

TRW has studied several digital correction techniques and evaluated 
these techniques in four primary areas: 

• processing-induced distortions 

• reseau and ground control point (GCP) location accuracy 

• Geometric distortion correction accuracy 

• throughput 

Processing-induced distortions include resolution degradation, corrected 
image discontinuities, and any artifacts. It is emphasized that this 
study was explicitly not concerned with image enhancement techniques; con- 
sequently, only resolution degradation and artifacts induced by the data 
processing were considered. The geometric distortion correction accuracy 
includes residual sensor distortions and attitude/ephemeris error-caused 
distortions. These are affected by the reseau and GCP location accuracies. 

Due to the shutdown of the RBV shortly after ERTS-1 launch, pro- 
cessing priority was shifted to MSS data although considerable progress 
had already been achieved in adapting TRW algorithms to accommodate the 
RBV data. As a consequence, the bulk of this report is focused on the 
evaluation of MSS data processing techniques. Table 1-1 lists the images 
processed in the course of the study. 
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Table 1-1. NASA Imagery Processed in the Course of the Study 


Sensor 

NASA 

Scene No. 

Geographical Description 

RBV 

1002-18134 

Monterey Bay 

MSS 

1057-18172 

San Pablo Bay 

MSS 

1062-15190 

Baltimore/Washington 

MSS 

1080rl5192 

Bal ti more/Washi ngton 
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1.2 STUDY APPROACH 

Evaluation of ERTS digital image correction techniques has proceeded 
in four stages: 

a. Adaptation of TRW algorithms to the ERTS data formats; 

b. Parametric design and performance evaluation; 

c. Performance evaluation of digital image corrections 
with data provided by NASA; and 

d. Determination of throughput for critical processing steps. 
Specific algorithm areas for evaluation included: 

a. Evaluation of three Interpolation techniques* 

• nearest neighbor interpolation 

t bilinear interpolation 
t TRW Cubic Convolution Process 

b. Kalman Filtering for precision attitude refinement 
(MSS only) 

★ 

c. Piecewise low-order distortion model 

d. Spiral search and shadow casting for reseau location * 

e. Automatic GCP extraction via sequential similarity detection 

In general, three approaches were employed for the evaluations: (a) auto- 

matic computer processing; (b) manual verification; and (c) visual inspection. 
Thus, evaluations of the geometric precision were assisted by computer pro- 
cessing for the rms differences between known ground locations and those 
determined from the corrected image (MSS and RBV), or by the rms differences 
between known reseau locations and those determined from the corrected 
image (RBV). 

Processing-induced distortions, which manifest themselves as resolu- 
tion degradation, image discontinuities and photometric errors were evaluated 
by means of statistical algorithms (histograms), spatial frequency analyses 
and visual inspection. Analyses of errors inherent in the image correction 
algorithm also will be presented. 

♦Developed prior to the commencement of this contract under company-sponsored 
programs . 
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Ground control point and reseau identification and location accuracy 
were evaluated by means of manual inspection and verification of photo- 
graphic and printer/plotter imagery for the areas in question. The speed 
and effectiveness of algorithms for GCP identification were also considered. 
Finally, the overall throughput of the processing was evaluated for RBV 
and MSS imagery. Comparisons were made between several resampling algor- 
ithms and examples of running times are given for CDC 6400 computer pro- 
cessing. 

The actual work progressed through three phases as illustrated in 
Figure 1-1. The first phase (pre-launch) was devoted to processing tool 
collection and modification. These tools were examined in the second 
phase with first-look ERTS-1 data, required tool modifications were initiated, 
and the compatibility of TRW's processing/evaluation techniques, with ERTS-1 
data was established. The second phase was concluded with the publication 
of the Data Analysis Plan for Phase III, the continuing data analysis phase. 
Phase III consisted of two principal facets: a) precision processing of 

selected imagery, and b) evaluation of the spatial and photometric proper- 
ties of the processed imagery from the standpoint of precision and through- 
put. Results of the Phase III studies are documented in this report. 

A top-level functional flow of the software tools utilized by TRW 
to process RBV and MSS data appears as the two step procedure illustrated 
in Figure 1-2. The first step, referred to herein as Pass I, includes 
the image distortion estimation function, the data reformatting function 
for bulk CCT data supplied by NASA/NDPF, and the GCP/reseau search region 
extraction function. The distortion estimation function results in the 
generation of a data file of distortion coefficients, corresponding to the 
piecewise bilinear distortion model employed. The second data pass corrects 

the reformatted bulk CCT data pixel by pixel, using one of three processes 
identified at the beginning of this section, making use of the distortion 
coefficients generated in the first pass processing. 

1 . 3 SUMMARY OF MAJOR RESULTS 

TRW has employed highly efficient digital techniques to precision 
process bulk ERTS data. These techniques were utilized for: 
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Figure 1-2. Top-Level Functional Flow for ERTS Data Processing 
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• Automatic GCP Extraction 

• Automatic Reseau Extraction (RBV only) 

• Kalman Filter Attitude Refinement (MSS only) 

• Piecewise Low Order Distortion Calculation 

• Image Correction. 

Both accuracy and throughput have been evaluated. 

Automatic GCP Extraction 

TRW has implemented an efficient technique which rapidly rejects 
mismatches between a reference Image subarea (32 pixels x 32 lines) and 
a test segment of the same size derived from a larger search region (150 
pixels x 120 lines). No positional errors were found when the reference 
subarea and search region were derived from the same ERTS bulk scene data. 
For search and reference data derived from two different orbital pass 
(i.e., 18 days or more apart in time), the maximum positional error was 
found to be 1 pixel. Running time on a CDC 6000 series computer averages 
about 14 sec. /GCP. 

Automatic Reseau Extraction 

TRW has implemented a reseau search technique which searches most 
probable reseau locations first. This is done in two ways: 1) ordering 

reseaux to be extracted in a sequence corresponding to increasing image 
distortion, i.e., center reseau in each row first, followed by next ad- 
jacent pair of reseaux, etc.; 2) spiral search of each reseau search region 

t 

(240 pixels x 240 lines) initiated at the position the reseau was last 
found (in the same or adjacent row), or at the position expected from pre- 
flight calibration measurements. This search is followed by thresholded 
shadow casting (row and column summing of image data) to locate the reseau 
center. Precision to a fraction of a pixel is derived from a polynomial fit 
to the shadow values. Incrementing the shadow threshold permits extraction 
of even low contrast reseaux. 

All reseaux extracted were correctly located, to the accuracy of 
independent measurement techniques utilizing line printer output. No false 

detections were experienced. Average time to extract a reseau was % 4 sec 
iof CPU on a CDC 6400 computer for a search region 240 pixels x 240 lines. 
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Kalman Filter Attitude Refinement 

The Kalman Filter was implemented in the form of an optimum linear 
attitude estimator, driven by errors in the UTM plane between actual GCP 
locations (obtained from the U.S. Geological Survey) and estimated GCP 
locations, using the BIAT data for Initial attitude estimates. Only 3 
GCPs were required to reduce residual errors to one pixel, a limiting value 
determined by the GCP extraction precision and sensor IFOV. CPU time for 
a CDC 6400 was 4.1 sec. 

Piecewise Low-Order Distortion Calculation 

For the MSS data it was found that a piecewise bilinear distortion 
model fit to 81 image positions (for which distortions were computed to 
high accuracy) is consistent with absolute worst case modeling errors <_ .5 
pixel. Highly accurate models for the sensor scan and spacecraft attitude 
and ephemeris are utilized to compute the resultant distortion coefficients 
for 8x8 regions spanning a scene, 8 coefficients per region. Distortion 
coefficient calculation for one MSS scene requires 20 sec of CPU exclusive 
of GCP extraction. Output corrected imagery is projected in a UTM coordinate 
system. 

RBV image distortions were also successfully modelled by means of a 
piecewise bilinear model, utilizing 36x36 search regions. Modeling errors 
in this case were <_ 1 pixel, and CPU time was about 44 sec, inclusive of 
reseau extraction, and exclusive of GCP extraction. 

Image Correction 

Precision correction of image data was accomplished by means of: 

1) nearest neighbor interpolation, 2) bilinear interpolation, and 3) TRW 
Cubic Convolution Process, developed by TRW prior to the commencement of 
this contract. It was found that both nearest neighbor and bilinear inter- 
polation processes significantly degraded image quality, whereas TRW Cubic 
Convolution generated the most satisfactory results. Thus, for example, 
nearest neighbor interpolation produced one pixel discontinuities, notice- 
able throughout the image. Bilinear interpolation was free of the image 
discontinuities produced by nearest neighbor interpolation, but noticeably 
degraded image resolution (spatial frequency content). Image data processed 
by the TRW Cubic Convolution Process had neither the nearest-neighbor dis- 
continuities nor the loss of spatial frequency content occasioned by bi- 
linear interpolation. Furthermore, comparison of TRW's digitally processed 
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ERTS Imagery with current NASA precision products reveals the superior 
quality of the digital processed result to the analogue method (scan/ 
rescan) . 

CPU time for the nearest-neighbor interpolation process was fastest 
(about 25 sec/10 6 pixels), was intermediate for bilinear interpolation 
(about 56 sec/10 6 pixels), and for TRW Cubic Convolution was slowest 
(about 110 sec/10 pixels) for the three processes considered. This 
corresponds to a range of 12.5 minutes - 55 minutes for a full ERTS 
scene (4 bands, or a total of 30 x 10 6 pixels of output data) using a 
CDC-6400. 

1.4 CONCLUSIONS AND RECOMMENDATIONS 

TRW has developed and demonstrated all-digital techniques to pre- 
cision process ERTS imagery. Extension of the large-scale computer 
implementations to small-scale "minicomputer" configurations should 
ultimately result in a highly cost effective implementation (in fact 
at least one minicomputer offers a factor of two improvement in through- 
put over the CDC 6400 system utilized in this study). It is recommended 
that this approach be pursued, and that simultaneously efforts be under- 
taken to extend the existing techniques to permit registration of data from 
different 18-day S/C passes to a fraction of a pixel. This will make possibl 
change detection studies directly from digitally processed data. 

1.5 PUBLISHED PAPERS 

References (5), (6), (7) cite the papers published during the 
course of this contract. 

1.6 NEW TECHNOLOGY 

No new technology was discovered during the course of this con- 


tract. 
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2 

MSS DATA 
PROCESSING 


2.1 INTRODUCTION 

The Multi spectral Scanner Subsystem (MSS) gathers data by imaging 
the surface of the earth in several spectral bands simultaneously through 
the same optical system. The MSS for ERTS-1 is a 4-band scanner operating 
in the solar-reflected spectral band region from 0.5 to 1.1 micrometer 
wavelength. It scans cross- track swaths of 185 kilometers (100 nm) width, 
imaging six scan lines across in each of the four spectral bands simul- 
taneously. The object plane is scanned by means of an oscillating flat 
mirror between the scene and the double-ref lector, telescope type of 
optical chain. The 11.56 degree cross-track field of view is scanned as 
the mirror oscillates + 2.89 degrees about its nominal position as shown in 
Figure 2.1 (Reference (2)). Table 2-1 lists the spectral bands for ERTS-1 . 

The instantaneous field of view (IFOV) of each detector subtends an 
earth-area square of 79 meters on a side from the nominal orbital altitude. 
Across scan separation of IFOV ' s is approximately the same, but the separa- 
tion along scan is approximately 56 nm. 

Table 2-1. ERTS-1 MSS Band Assignment 


Band 


Wavelength 

Band 

1 

0.5 to 0.6 micrometers 

Band 

2 

0.6 to 0.7 micrometers 

Band 

3 

0.7 to 0.8 micrometers 

Band 

4 

0.8 to 1.1 micrometers 
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The sampling Interval is 9.95 ysec, corresponding to a cross -track 
motion of the instantaneous field of view of 56 meters. The along-track 
scan is produced by the orbital motion of the spacecraft. By oscillating 
the mirror at a rate of 13.62 Hz, the subsatellite point will have moved 
474 meters along-track during the 73.42 millisecond active scan and re- 
trace cycle. The width of the along-track field of view of six detectors 
is also 474 meters. Figure 2>1 shows the composite scan pattern. 

The downlinked MSS data is processed by the NASA Data Processing 
Facility (NDPF) to produce film products and Computer Compatible Tapes 
(CCT). Four CCT's are required for the digital data corresponding to one 
scene, framed 185 Km on a side. The MSS data on the CCT's is spectrally 
interleaved, line-length adjusted, and radiometrically calibrated. For a 
given spectral band, one active scan of the MSS mirror produces simultan- 
eously six lines of bulk image data. A total of 2340 lines (390 scans) 
constitute one such image. For each set of images there is a set of S/C 
ephemeris and attitude data, relative to the time of the image center 
point (image center time, ICT), recorded on the bulk image annotation 
tape (BIAT) . 

TRW's first pass program for multi spectral scanner data is called 
MSS. MSS is a CDC 6400 program encoded in modular manner to flexibly 
exercise PASS-I algorithms. 

2.2 METHODOLOGY 

The function of the MSS program is to generate a file of distortion 
coefficients representing the geometric correction transformation which 
maps the bulk image to the geometrically corrected output image in a UTM 
projection. The geometric distortions of the MSS bulk image are due to 
errors in the vehicle attitude/ephemeris which may vary from scan to 
scan, to irregularities in the motion of the scanning mirror, and 
to earth curvature and rotation, A schematic of the MSS software is 
given in Figure 2.2. 

The determination of the distortion coefficients requires accurate 
models for the sensor scan and the vehicle attitude/ephemeris. Initial 
values of spacecraft attitude, altitude, and nadir attitude and longitude 
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Figure 2.2 MSS Program Functional Schematic 
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are provided on the annotation tape for 9 equally-spaced points spanning 
an image. These values are least-squares fit to corresponding cubic poly- 
nomials in time to develop an initial estimate of the vehicle attitude/ 
ephemeris. The attitude data is then refined, on option, by determining 
positions of selected Ground Control Points (GCP's) on the input CCT, and 
employing a Kalman filter for a precision refinement of the roll, pitch, 
and yaw coefficients. 

Given the refined (or initial) attitude profile, the image distortion 
coefficients are computed by first defining a regular grid of image loca- 
tions, called "pseudo-reseaux," in the coordinate system of the bulk image. 
These points are mapped to their corresponding locations in UTM coordinates 
in the corrected image coordinate system. The LOS vector from the space- 
craft is determined for each pseudo reseau point as a function of time 
through the use of analytical models of the earth-vehicle system and the 
MSS mirror scan. 

A regular grid of pseudo-reseaux is then defined in the UTM system 
parallel to the spacecraft ground track and centered with respect to the 
initial set of pseudo-reseaux mapped from the bulk image system by means 
of a least-squares fit of the edge reseaux. The distortion coefficients 
(defined for each block of four pseudo-reseaux in the regular UTM system) 
are computed on the basis of a bilinear distortion model, through an in- 
verse mapping from the regular reseau coordinate system to the bulk image 
(uncorrected) coordinate system. This mapping is accomplished by means of 
a rapidly convergent iterative technique. 

2.3 MSS DATA PROCESSING RESULTS 

The results of the evaluation of TRW algorithms for processing MSS 
data encompass five areas: GCP extraction precision, Kalman filter per- 

formance, Pass I precision, resampling algorithm comparisons, and 
module throughput. 

2.3.1 GCP Extraction Precision 


GCP search regions called "neighborhoods" (see Figure 2.3) 


were 
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Pixels 


Lines 



Figure 2.3. GCP Window Extraction from a Search Neighborhood 

Dimensions refer to line and pixel values in bulk CCT data. A 
window is located within a neighborhood when the cumulated 
pixel value differences (random pixel ordering) within an 
overlap region (shown dotted above) do not exceed a moving 
threshold. A possible location is rejected quickly when the 
cumulated difference exceeds the same threshold, after a pre- 
determined number of differences are cumulated. 
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extracted as part of the bulk data reformatting process*, which produces 
line sequential and band separated bulk image data. GCP location within 
each search region is accomplished by means of a feature matching technique 
(Reference (2)) which rapidly rejects mismatches. 

Typical of the precision of the technique are the results given in 
Table 2-2. Four tests were made using scene 1062-15190-5 and scene data 
obtained 18 days later, 1080^15192^5. For simplicity, let scene 
1062-15190-5 be denoted as 62, and let scene 1080-15192-5 be denoted as 
80; let N (neighborhood) represent the search retion, (150 pixels x 
120 lines), and let W (window) represent the reference template of the 
GCP (32 pixels x 32 lines). 

The intersection of Interstate 83 and the Baltimore Beltway stands 
out clearly in both scenes and was taken as the GCP. The tests consisted 
of utilizing: a) 62 as the window data (62W), and 80 for the source of 

neighborhood data (80N); b) 62 as both sources for window and neighbor- 
hood data (62W + 62N); c) 80 as both sources for window and neighborhood 
data (80W + 80N); and d) the reverse of the first case (80W + 62N). The 
paramater XN0ISE is the threshold increment used to compute the threshold 
0T. Basically, a large XN0ISE will give a large probability of 
false detection and cause the algorithm to run slower, while a small 
XN0ISE will give a large probability of false rejection and cause 
the algorithm to run faster. For XN0ISE=4 , the mean time to extract a 
GCP is 19.2 cpu seconds, and the errors for 62W + 80N and 80W + 62N 
reach a minimum value. 

Observations made subsequent to these tests indicate that a one 
pixel error can occur if the feature in question has been nearest-neighbor 
line length adjusted, i.e., a replicated pixel is present (Reference (3), 
Appendix 2). 


Using CDC computers necessitated conversion to 6 bit data per 
pixel at the same time. 

Out of a possible range of 64, for 6 bit data. 
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Table 2-2. Typical GCP Extraction Precision 


WINDOW + 
NEIGHBORHOOD 

XNOISE 

TIME TO 
EXTRACT GCP 
(cpu sec) 

ERROR 

LINES PIXELS 


1 

4.8 

0 

0 


2 

5.5 

0 

0 

62W + 62N 

3 

7.9 

0 

0 


4 

9.6 

0 

0 


1 

5.1 

0 

0 


2 

7.2 

0 

0 

80W + 80W 

3 

11.2 

0 

0 


4 

25.0 

0 

0 


3 

6.8 

18 

32 


3.5 

8.3 

-28 

-36 

62W + 80N 

4 

13.3 

-1 

-1 


5 

44.1 

-1 

-1 


6 

> 60 

-- 

-- 


3 

8.9 

12 

36 


4 

15.0 

0 

0 

80W + 62N 

4.5 

34.5 

0 

0 


5 

> 60 




NOTE: In no case was an error in GCP location encountered when the 

GCP window was extracted from the neighborhood in question, 
Error sources for neighborhoods and windows derived from 
different scenes include: a) scene photometry differences, 

related to the seasonal variation in solar elevation angle; 

b) changes in vegetation with season (brown and green waves); 

c) line length correction differences between two scenes, in- 
volving the nearest-neighborhood replicated pixel occurring 
in the feature for only one of the two scenes. 
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2.3.2 Kalman Filter Performance 

The precision rectification of an ERTS/MSS bulk image requires an 
estimate of the ERTS attitude time-series over the time interval in which 
the image was scanned. Such an estimate can be obtained by fitting poly- 
nomials to the sampled attitude data on the Bulk Image Annotation Tape 
(BIAT). Unfortunately, an attitude time-series obtained in this way does 
not meet the stringent requirements of precision image rectification. To 
achieve image rectification on the order of 1 pixel, it is necessary to 
know each attitude component to within +.01 milliradians. According to 
the ERTS Data User's Handbook (Reference (1)), the ERTS attitude measure- 
ment system is accurate only to the nearest milliradian. Clearly, some 
additional information is required to obtain a more refined estimate. 

Since the attitude data available from the annotation tape cor- 
responds to a limited number of values, the MSS program utilizes a Kalman 
filter to update the initial estimates of attitude coefficients using 
Ground Control Points. When this option is selected, the Kalman filter 
subroutine updates the attitude coefficients given the initial coef- 
ficients (state vector), the error covariance matrices for each component, 
and the differences in GCP ground locations (known a priori location minus 
that computed from the bulk image). 

The estimator was implemented on the CDC 6400 computer using data 
derived from ERTS/MSS bulk image, 1062-15190 ( 1 068-35 : BIAT) . For the image 
considered, it was shown that only 3 Ground Control Points are required 
to drive the Ground Control Point position estimation error from 4.5 
kilometers to less than 80 meters. This result suggests that the esti- 
mator is performing up to the resol utional capability of the MSS detectors. 
GCP 's were selected from a list supplied by USGS and 7.5 minute quadrangle 
maps provided by USGS. References (6) & (8) describe TRW's 2-dimensional 
UTM observable estimator. Let the ERTS attitude time-series (restricted to 

the image time interval) be represented by e = (e , e , e ) whose sub- 

^ -p 

vector components are comprised of sets of coefficients which define 
polynomial realizations of the roll, pitch, and yaw time series. Let a 
represent how well the ground truth points are known. 
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The prior estimate of e is derived from attitude-sampled-data on 
the BIAT. Figures 2.4, 2.5, and 2.6 show the results of least-squares 
fitting 1 st , 2 nd , and 3 rd order polynomials to each set of sampled -data. 
These figures demonstrate the sufficiency of a third-order polynomial 
representation for each of the attitude components. 

Due to the unavailability of statistics, the a priori covariance 
matrix P(o) and the GCP location accuracies o(k), k=l,2,..., were assigned 
on a subjective basis. The selection of P(o) was circumscribed by the 
following assumptions: 

1) The components of are uncorrelated. 

2) The uncertainty in the zeroth order component of e 
(ae{r,p,y>) is bounded by the attitude measurement 01 
uncertainty specified by NASA (1.2217 milliradians) . 

3) Uncertainty in the j tl1 order component of e 
(j=l,2,3; ae{r,p,y}) does not exceed 10% _a 

of the absolute value of the component's initial value. 

On the basis of these assumptions, P(o) was taken to be diagonal with 

( P ( 0 ) ) i i =S (i^l )mod(4)+l 1=1,2 ****» 12 

S 1 = 1 . 221 73x1 0~ 3 , S 2 = 10~ 7 , S 3 = 10" 8 , 

S 4 = 10“ 10 . 

For GCP's not obscured by cloud cover, o was assigned a value of 30 
meters. 

The filter representation shown in the appendix is that of an 
operational implementation. Classically, the state update is written: 

*NEVT *0LD + " ^TRUE^ 

where A is the gain matrix and the ^ vectors are GCP locations. The norm 
of the GCP residual error, y. - ^-true* was founcl t0 a usef ul measure of 
filter stability. Intuitively, the norm should decrease with increasing 

a v 

observations of GCP's if e_ (estimate of e) is stable. Figure 2.6 indicates 
this behavior for 10 GCP's. 


21 7'L- 
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Figure 2.4. Polynomial Models of Roll Time-Series 



2nd-order 
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Figure 2.5. Polynomial Models of Pitch Time-Seri 


2nd-order 
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Figure 2.6, Polynomial Models of Yaw Timer-Series 


Initial and Final Estimates of Attitude Polynomial Coefficients 
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Figure 2.7. Norm of GCP Residual Error Vs GCP Observation Number 



20634-6003-TU-02 
Page 2-15 


th 

Another measure of e_'s stability is the behavior of the zero 
order error-covariance terms as a function of observation number. Plots 
of the square roots of these functions (Figure 2.8), exhibit a highly 
stable monotone decreasing behavior which tends to be asymptotic to zero. 
Notice that after 10 GCP's, the roll and pitch estimates have converged to 
within +_ .01 miniradians of their true values, whereas the yaw estimate 
converges to within + .28 milliradians of its true value. 

Comparison of the norms before and after filtering provide a measure 
of the filter's worth. Before filtering there is a persistant error on 
the order of 4.5 kilometers; but after filtering, the error is on the 
order of 79 meters (1 pixel), suggesting that the estimator is performing 
to the limits of the sensor IFOV and GCP location accuracy. 

Figure 2.9 shows the accuracy of the estimator as a function of the 

number of observations, k=l , 2 10. After the k th observation, the 

resulting e_(k) was used to compute the norm for all observations and a 
mean standard deviation were computed (sample size=10). The plots in 
Figure 2.8 show the results of repeating this process for k=l , 2, ...» 10. 
The figure shows that only three GCP's are required to correct the 
attitude time-series to one GCP, since the attitude uncertainty is 
zero*”* 1 order, as evidenced by comparing the initial and final estimates 
of £ given in Figure 2.7. 

2.3.3 Pseudo-Reseau Grid Modeling Errors 

Corrections for geometric distortions in the MSS image are obtained 
by calculating exact corrections for a subset of image points (pseudo- 
reseau points) and applying a piecewise bilinear distortion model to 
determine corrections for all other points. The accuracy of the piece- 
wise model increases as the number of pseudo-reseaux is increased; how- 
ever, processing time also increases. The relationship between number of 
pseudo-reseaux and positional error was investigated to permit selection 
of pseudo-reseau grid size to meet a desired error criterion. 
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The MSS Pass I program was modified as follows: 

1) The Ground Control Point (GCP) and Kalman filter 
routines, ordinarily used for attitude refinement, 
were deleted. 

2) A provision for inputting "test points" was added 
(a test point is a pixel for which the error is to 
be determined). 

3) The exact projection of each test point onto the output 
space was determined using the existing subroutines. 

4) An additional subroutine was added to perform the bilinear 
interpolation for arbitrary pixel locations. 

5) The difference between the exact projection and inter- 
polated position of each test point was taken as the 
error. 

6) Summary statistics for the errors are calculated. 

Modeling errors were studied by (1) systematically scanning the image 
to locate the actual maximum error; and (2) randomly selecting test points 
(N=200) to obtain frequency distributions of errors as well as values for 
mean errors and standard deviations. Several measures were calculated: 
mean error, mean absolute error, maximum error from the systematic scans, 
maximum error from the random selections, error standard deviation a, 
mean absolute error +2 a, and the error frequency distribution. Errors 
were measured separately in the line and pixel directions. 

The principal results are shown in Table 2-3 for pseudo-reseau grid 
sizes ranging from 5 x 5 to 13 x 13 points. If a maximum (worst case) 
allowable position errors are on the order of 0.5 pixel or 0.5 line, the 
results indicate that a 9 x 9 pseudo-reseau grid is required. Further, 
it seems unlikely that a grid smaller than 7x7 (maximum errors on the 
order of 1.0 pixel or line) or greater than 11 x 11 (100% of errors less 
than 0.4 pixel or line) would be considered. Note that mean absolute 
errors are substantially smaller than the maximum errors and that the 
absolute error mean +2 a values (+ 2 a about the mean is a 95% confidence 
limit, assuming a normal error distribution) agree fairly well with the 
actual maximum values obtained, whether from the random selection or 
systematic scan methods. 
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Table 2-3. Position Error Data for Scene 1057^18172 


PSEUDO-RESEAU 
GRID SIZE* 
(AS) x (SA) 

TYPE 

r~wm 

ABSOLUTE 

ERROR 

o 

MAXIMUM ERROR 
(RANDOM/SYSTEMATIC) 

MEAN 1 

+ 1 

2a 

5x5 

PX 


.49 

1.26 

1.38 


LN 

.43 

.55 

1.28 

1.53 

7x5 

PX 

.21 

.26 

0.64 

.85 


LN 

.41 

.51 

1.26 

1.43 

7 x 7 

PX 

.23 

.26 

.54 

.75 


LN 

.22 

.28 

.73 

.78 

7x9 

PX 

.21 

.26 

.79 

.73 


LN 

.11 

.15 

.39 

.50 

9x7 

PX 

.13 

.13 ! 

.39 

.39 

j 

LN 

i 

i 

i 

.22 

.28 

.71 

.78 

9x9 

H3|H 

MM 

ill 

.55 

■Ml 



Mm 

wm 

.39 

■Ml 

11 x 9 

PX 

.11 

.14 

.47 

.39 

! 

! 

LN 

.11 

.14 

.39 

.39 

s ! 

11 X 11 

PX 

wm 

.11 

1 

.33 

.36 


LN 

1 

.10 

.38 

.29 

13 x 13 

PX 

.08 

.10 

.29 

.28 

i 

LN 

.06 

.08 

.26 

.22 

i 


NOTE: Grid dimensions are given as Along Scan (AS) x Scan Advance (SA). 

Errors are given in the Along Scan (PX) and Scan Advance (LN) 
directions. 
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(The a value was obtained from the algebraic error values.) The larger 
maximum random error or maximum systematic error is shown in Table 2-3. 

A plot of the mean absolute errors and maximum errors as determined 
from the random test point runs is shown in Figure 2.10 for scene 1057-18172. 
The actual number of reseau blocks In the scan advance and along scan 
directions is given on the abscissa. In some cases, there are several 
points for a given number of reseau blocks as a given number of blocks 
in one direction was repeated while the number in the other direction 
was varied (e.g., 7x7 and 7x9 cases from Table 2-3. The along scan 
points show more variability than do the scan advance points, indicating 
that along scan interpolation errors for a given number of along scan 
blocks might be more sensitive to the number of scan advance blocks than 
vice versa. These curves can be used to predict (approximately), the 
error values for reseau sizes not included in this study. 

Three grid sizes were run for a second picture, scene 1057-16375. 

The results indicate that the pixel error is about the same as for 
scene 1057-18172, but that the line error is greatly reduced. In general, 
maximum errors were found at the image corners. Within a pseudo-reseau 
block, maximum errors usually occurred in the vicinity of the center of 
the block. 

Other variations were tried, including an azimuthal rotation of 
0.5°, and an offset in the input grid origin from which the pseudo-reseau 
points are specified. None had an effect on error values. In summary, 
acceptable errors can be obtained with a reasonable number of pseudo- 
reseaux. The selection of a reseau grid will depend on the error size 
that is acceptable for a given purpose; a 9 x 9 grid would seem to be a 
good choice based on the images studied. 


ADVANCE (LINES) 



Figure 2. IQ. Mean Absolute and Maximum Random Position Errors 
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2.3.4 Pass I Precision 

Using the results of the previous sections, it is possible to 
establish an error budget for the Pass I process. This budget is sum- 
marized in Table 2-4. Error sources modeled by the Pass I software ele- 
ments include: spacecraft altitude variations within a frame; spacecraft 

attitude variations (roll, pitch, yaw) within a frame, earth rotation; 
spacecraft velocity changes; finite scan time effect; scan mirror velo- 
city nonlinearities; perspective corrections; and bias errors due to 
MSS alignment to AMS. 

Note that the 6CP extraction error, while listed in the table, is 
not included in the computation of the RSS/RMS error. It is excluded 
since the GCP extraction error directly relates to the Kalman filter 
error-- the filter will do no better than the GCP extraction accuracy. As 
a consequence, the Pass I RSS error is 1.12 pixels and the Pass I RMS 
error is 0.79 pixels. 


Table 2-4. Pass I Error Budget 


PASS I ELEMENT 

MAXIMUM ABSOLUTE ERROR 
(Pixels) 

GCP Extraction 

1.0 

Kalman Filter 

1.0 

Pseudo-Reseau Modeling 

' 

0.5 

RSS 

1.12 

RMS 

0.79 
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2.3.5 Processed Imagery 

Typical of MSS image processing performance is scene 1062-15190, con- 
taining tne cities of Baltimore and Washington, D.C. Figure 2-11 shows the 
bulk image for band 5 (red), the original negative for which was reproduced 
directly trom NASA-supplied bulk CCT data, with no correction tor aspect. 
Figure 2-12 is a NASA precision processed image for the same scene and band.* 
The loss of information proceeding from bulk CCT data to the NASA product is 
quite evident, and shall be made even clearer by the following figures. 

Figures 2.13 and 2.14 show a detail (for bands 5 and 7, respectively) 
resulting from TRW's precision digital processing techniques and TRW's Cubic 
Convolution Process. The detail (scale 'v 1:571,000**) is centered around 
the city of Baltimore and its harbor. This region was selected instead of 
Washington because of the richness ot detail and because silt in the Poto- 
mac resulting from storm activity reduced land/water contrast significantly 
(in the visible). 

It is evident from Figure 2.14 at once that there is a considerable 
amount of detail in the digitally processed product. Note for example the 
runways of a Baltimore-Washington International Airport (Friendship Inter- 
national), which are not evident in the NASA-precision product, an enlarge- 
ment of which (band 7) is shown in Figure z.15. 

For comparison purposes. Figure 2.16 shows bulk image subareas corres- 
ponding to Figures 2.13 and 2.14 (bands 5 and 7), reproduced trom bulk CCT 
data. Also for comparison purposes, Figure 2.17 shows a full scene precision 
product generated by TRW's digital methods. It can be compared with the 
NASA product shown in Figure 2.12. 

Additional MSS imagery was previously processed and reported on in the 
first 6-month report for this contract (Reference (8)) and in Reference (5). 
Comparisons between and analyses of three resampling methods are contained 
in Section 4 of this report. References (8) and (9). 


The original for this image was created as a contact negative from the 
NASA-supplied positive 9 1/2 inch transparency. 

On the original negative 100p = 57.098765m. 
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Figure 2.11. Bulk Image 1062-15190-5 from CCT Data 
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Figure 2.12. NASA Precision Product for Image 1062-15190-5 
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Figure 2.13. Detail of TRW Precision Product for Image 1062-15190-5 
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Figure 2.17. TRW Precision Product for Image 1062-15190-5. 

This product was generated on a precision filmwriter from 
digitally processed bulk CCT data supplied Dy NASA. 
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2.3.6 Throughput 

Table 2-5 summarizes the throughput characteristics for the several 
Pass-I MSS modules and the Pass-II algorithms implemented in FORTRAN on 
TRW's 6400. Times are given as central processor seconds. Pass-I executes 
under 17 seconds excluding GCP refinement, and less than 1 minute with 
GCP extraction and refinement (3 GCP's). The Pass-II software, operating 
on an image of 10^ pixels, requires 52 seconds for bilinear interpolation 
and 110 seconds if the TRW Cubic Convolution Process is selected. The 
Pass-II nearest-neighbor interpolation requires about 30 seconds/10^ 
pixels. 

Cursory extrapolation of these numbers to a minicomputer indicate 
that Pass-I should execute as fast or faster than the CDC FORTRAN 
implementation. The bulk of the computing load, however, is in the 
Pass-II process; minicomputer throughput for the bilinear software is 
estimated at less than 5 minutes/MSS image, while the TRW Cubic Convolu- 
tion Process is estimated to require less than 10 minutes/MSS image. 

Table 2-5. Module Throughput 


ELEMENT 

CPU TIME (sec) 

PASS-I 


1. Attitude/ephemeris time- 
series fittings; distortion 
coefficient calculation 

16.1 

2. GCP exptraction )per GCP 

14.2 

3. Kalman filter refinement 
(15 GCP) 

4.1 

PASS-II (10 6 Pixels) 


1. Nearest-neighbor 

30 

2. Bilinear 

52 

3. TRW Cubic Convolution 

no 
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3 

R B V DATA 
PROCESSING 


3.1 INTRODUCTION 

The Return Beam Vidicon (RBV) camera subsystem contains three 
individual cameras that operate in different nominal spectral bands. The 
measured spectral response of the three cameras is shown in Figure 3.1. 

Each camera contains an optical lens, a shutter, the RBV sensor, a thermo- 
electric cooler, deflection and focus coils, erase lamps and the sensor 
electronics. The cameras are similar except for the spectral filters con- 
tained in the lens assemblies to provide separate spectral viewing regions. 
The sensor electronics contain the logic circuits to program and coordinate 
the operation of the three cameras as a complete integrated system and 
provide the interface with the other spacecraft subsystems. 

The three RBV cameras are aligned in the spacecraft to view the same 
nominal 185 kilometers (100 nautical mile) square ground scene as depicted 
in Figure 3.1. When the cameras are shuttered, the images are stored on the 
RBV photosensitive surfaces, then scanned to produce video outputs. The 
three cameras are scanned in sequence during the last 10.5 seconds of the 
basic 25 second picture time cycle. The video from each is serially com- 
bined with injected horizontal and vertical sync. The readout sequence is 
camera 3, then camera 2, then camera 1. 

The digital reduction of the RBV data involves the following proces- 
sing: 

• Search for and measure positions of reseaux 

• Perform a polynomial fit to measured reseau positions 
and extrapolate positions of lost reseaux 

• Search for ground control points and measure displace- 
ments from nominal positions 
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RETURN BEAM 
VIDICON CAMERAS (3) 




SOURCE: REFERENCE 2 


Figure 3.1. RBV Subsystem 





20634-6U03-TU-02 
Page 3-3 


• Compute grid of pseudo-reseaux positions at 113 pixel 
intervals based on biquadratic coefficients determined 
from displacements of 9 nearby reseaux, and ground con- 
trol point displacements 

• Compute piecewise bilinear distortion ceofficients for 
intervals between pseudo-reseaux 

• Resample entire scene on a regular output grid (see Section 
2.3.5 for discussion of resampling techniques used in this 
study) . 

Because of the large amount of digital data in each scene, those 
processes involving comparisons or interpolation at individual pixel levels 
have been analyzed and coded for high efficiency. The search for the 
reseaux, the search for the ground control points, and the final resampling 
require pixel by pixel treatment of large segments of the data. As a re- 
sult, sophisticated analytic techniques have been used to reduce the number 
of pixels handled (i.e., the searching is made efficient by using as much 
information and insight as possible). 

3.2 APPROACH 

3.2.1 Reseau Extraction 

The reseau marks for an RBV frame are placed in a 9 x 9 array on 
the RBV tube so as to be superimposed on the RBV image. They are placed 
with an accuracy of 3 micrometers or 1/4 pixel accuracy. The RBV face- 
place reseaux have three possible locations: (1) intended, (2) nominal, 

(3) measured (found). The intended locations of reseaux lie on a regular 
grid and represent the desired positions. Due to imperfections in the 
process by which they are created on the faceplate, their actual locations 
are slightly different from the intended locations. Calibration data 
supplied by NASA corresponds to the nominal reseau positions. The loca- 
tions actually found as a result of Pass 1 processing are denoted as 
measured positions. In addition, reference locations in the output cor- 
rected image correspond to regular grid locations (very close to the in- 
tended grid positions) with respect to which the Pass 2 coefficients are 
def i ned . 

The nominal displacement in each coordinate is less than 1 pixel; 
the measured displacement is normally less than 3% or 120 pixels. The 
dimensions of the entire RBV frame is 4125 lines with 4608 pixels per 
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line in bulk format. 

For each reseau a 240 x 240 pixel neighborhood is called from disk 
storage, with the intended reseau location at the center of the neighbor- 
hood. In searching for the reseau, a neighborhood is divided into blocks 
of 15 x 15 pixels. Excluding a band 15 pixels wide at the edges of the 
neighborhood, each block is searched to see if it contains the center of 
the reseau. Because a reseau is 32 pixels wide, 32 pixels high, and 
4 pixels across each arm, the center must be at least 16 pixels from the 
edge of the neighborhood for the whole reseau to be contained in the 
neighborhood. Therefore, the center cannot be in the outer 15 pixels, and 
only the inner blocks of the neighborhood need be searched for the reseau 
center. 

The order in which the blocks are searched is determined first by 
the block in which the last reseau was found. If this block and those im- 
mediately adjacent to it fail to contain the reseau center, then a spiral 
pattern of blocks is searched, moving outward from the intended reseau 
location (skipping those searched in the first try near the block pre- 
dicted from the previous RBV image). Before a block is searched for the 
reseau center, the border of the block is examined to see if any edge 
of the block fails to have a black pixel (checking every other pixel). If 
the center of the reseau is inside the block, then at least two adjacent 
pixels on each edge must be black because one of the reseau arms will cross 
it. The block is not searched if it cannot meet this test. 

To search a block for the reseau center, a check is made of every 
second pixel in alternate rows from the top to the bottom of the block. 

When a pixel is found to be black, then pixels at locations + 12 and + 6 
in both coordinates are also checked to see if they are black. (Note that 
only the center pixel is required to be in the block being searched. Some 
of the surrounding pixels tested will be outside the block, but still in- 
side the neighborhood. Hence, their intensities will be available.) 

After a candidate center has passed the search criteria described, a 
shadow casting routine is used to determine whether the real reseau center 
is within + 2 pixels of the candidate center. Only a 27 x 27 pixel reseau 
is sought (that is, the sums for the shadow casting include 27 pixels). 
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since the outer pixels may be lost In distortions, and so that the 
reseau will be detected even if the candidate center is displaced as 
much as two pixels from the true center. A set of 13 line sums is 
made in each coordinate, centered about the candidate center. Seven 
sums are considered in the search for a minimum intensity, while the 
outer 6 sums are available for information on the background. If one of 
the sums in each coordinate is below an input threshold (27 x IT) and if 
all of the sums three pixels away from it are at least an input increment 
(IDEL) greater, the reseau is found. The threshold increment ITINC (=7) 
for IT (=3) is cumulated until the sum of IT+ITINC's exceeds ITOP if a 
reseau does not pass this test. No more than 3 iterations (IT0P=17) were 
required in this study, and no false alarms were encountered. 

The precise center location is then found by a three point interpola 
tion around the minimum sum in each coordinate. Finally, the nominal 
displacement is subtracted from the measured displacement to obtain 
the distortion correction to be applied at the intended reseau location. 
The known reseau displacements are used to find the coefficients a, 
b, for the global third order polynomial: 

3 3-i . . . -1 

<Sx = E E a i . xy 1 a = (A T A) A T 6x 

i=0 j=0 

<$y = £ if b u xV b = (A T A) 1 A T sy 

1-0 j=0 


All reseau displacements are then predicted. The known reseau 
displacements are tested to see whether they differ an unacceptable 
amount (TH) from their predicted displacements. Those rejected, and 
those not found originally are assigned the predicted displacement value. 

An option permits the recomputation of the predicted displacements 
with the rejected reseaux eliminated from the coefficient determination 
(an iteration). In this case, any additional reseaux that fail to be 
close enough to the recomputed predicted value are simply assigned their 
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predicted value. 


3.2.2 Ground Truth Computation for Attitude and Ephemeris Error 

Because of vehicle attitude and position uncertainties, ground 
control points are located. Their displacements (6x^ , 6y.j, 1=1 . N) on 
the RBV image, are used to solve for the corrections to vehicle position 
and attitude: 


B 





o 


e 

e 

e 

6 


R 

P 

Y 

H 


where 6 x q , 6y Q , 6H are position errors along-track, cross-track, and in 
attitude, and e^, 0p, 0y are departures from nadir in the roll, pitch, 
yaw coordinates. Inputes estimates B Q for B, and its inverse covariance 
matrix, SI, obtained from the ephemeris tape data, are combined with the 
solution from ground control displacements to yield a best estimate, B 
of B: 


B = (M T M + si)' 1 (M T aW + SI * B q ), 


6W = 


6X, 


6X W 
4 *1 


m 


where 
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and where the rows of M are: 

M-j = C 1 * 0 , (Hx xH, jjp] 

for 1 _< i <_ n, and for n+1 <_ j <_ 2n 

Mj = [0, 1, (Hx xH, ^-] 

The search for the ground control points requires a fast and simple 
means to recognize known ground features. The method of Reference (3) 
has been coded, with modification, to allow the inclusion of a set of 
priority pixel points for comparison before the comparison of random 
points. A window of dimension 32 x 32 pixels about the desired ground 
control point, is assumed for the given site. As many as 50 pixels 
within the window may be entered as special priority pixels for scene 
recognition. A neighborhood of 120 x 120 pixels off the RBV image is 
taken from disk storage to compare with the window. By knowing the 
reseau displacement for a nearby reseau, the displacement of the ground 
control point on the RBV image from its nominal position due to distortions 
can be predicted, and the residue of attitude displacement and ephemeris 
error can be limited to +60 pixels). A larger block of data is stored 
in the initial stripping of data off the bulk tapes (300 x 300 pixels), 
and the predicted displacement used to enter the larger block to obtain 
the 120 x 120 pixel neighborhood. 

In brief, the technique used to recognize the ground control point 
involves repeated comparisons of the window with subsets of the neighbor- 
hood. Each comparison is made until the summed residuals between the 
compared pixels in the window and the subset exceed a threshold. The 
number of pixels compared before reaching the threshold is used as a 
measure of the goodness of fit between the two scenes. The threshold 
is an increasing one, based on the expected residuals due to noise. Both 
the window intensities and the subset intensities are normalized before 
differencing to obtain the residuals. The order in which pixels are com- 
pared is as follows: First the priority pixels are compared, and then a 

random order of points established one for all subset comparisons is 
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followed. The randomness assures that a variety of information is used 
to measure the fit. 

The summation of normalized residuals is as follows: 

Xft- E < s i - s - w, - W) 

where is the intensity for pixel i in the subset, is the intensity 
for pixel i in the window, S is the mean intensity for the subset, and W 
is the mean intensity for the window. Subroutine NORMAL computes A and W. 
Subroutine RAND assigns the order of i. Subroutine SURFAC performs the 
comparisons and keeps track of the subset having the largest number of 
comparisons before exceeding the threshold. The threshold is based on 
the results of Reference (2 ) and may be written as: 

JT = 9 x NOISE + NOISE 
i 

and is easily computed by a single add for each comparison. NOISE is the 
mean of the image noise. 

The normalization process requires that S be recomputed for each 
subset. This is accomplished by incrementing the previous value to 
avoid (32) 2 additions. 


3.2.3 Pseudo-^Reseau Grid Generation 

In a study performed by TRW prior to this contract, modeling error 
analyses were performed for piecewise biquadratic and piecewise bilinear 
distortion models, and global third order and biquadratic interpolation was 
found to be most efficient and accurate for the RBV problem. For the error 
resulting from the interpolation statistics themselves, due to spatial ly- 
uncorrelated distortions and reseau measurement errors, the piecewise 
biquadratic method results in a ratio of error in interpolation to error 
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In reseau position, a/a R , of 1 for points internal to the grid and 5 at 
the corners. Figure 3.2 shows graphically how the ratio varies over a 
scene. For example, if the error in measuring reseaux is 1 pixel, the 
internal error for the interpolation would be 1 pixel, and the error at 
the edges and corners of the scene would reach 5 pixels, resulting from 
the piecewise biquadratic distortion model used to obtain pseudo-reseau 
positions. (This biquadratic model is the major error source). The 
sensitivity of the piecewise biquadratic to modeling errors for 1% first- 
order pincushion, 0.2% second-order pincushion, and 0.5 pixel 1 -sigma 
random distortion (a bad set of distortions) will result in errors of: 

Up to 1 pixel at corners of scene 
Up to 1/6 pixel internal to scene 
1/4 pixel RMS error over all of the scene 
The other methods tested produced errors as much as 10 times larger. 

To utilize bilinear interpolation of distortions to compute the geo- 
centric correction for each pixel, a pseudo-reseau grid is computed using a 
global biquadratic polynomial. The polynomial predicts the displacement 
for pseudo-reseaux placed every 113 pixels between the reseaux, which are 
at 452 pixel intervals. At the outer edges of the reseau grid, two rows 
of pseudo-reseaux are extrapolated. Nine reseau points are used at a 
time to obtain the pseudo-reseaux contained between them and on the outer 
borders. Figure 3.3 shows the 16 reseau subregions of 9 reseaux each 
used to compute pseudo-reseaux. The original reseau points are included 
in the 37 x 37 matrix of pseudo-reseau points prepared for the final bi- 
linear interpolation. 

The coefficients for the global biquadratic fit are contained in: 

t •u-'V-iA 

i=0 i=0 

and are computed from 

-1 

a = A 6x 
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Figure 3.2. The Statistical Function (a/a R ) for the Piecewise 

Biquadratic Approach 
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Figure 3.3. Reseau Subregions 


Intended locations of reseaux define border of regions, with 9 reseaux 
in or on the border of each region for use in biquadratic interpolation 
within the region. The measured reseau displacement minus the nominal 
displacement is assumed to apply at the intended reseau location. 

For the ij th reseau, the intended location is 

M x = 452i -2260 
M y = 452J-2260 
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where the A matrix has rows 

P1 2 2 2 2 2 2 

[l»x,x ,y,xy,x y,y ,xy ,x y ] 

one row for each of the nine reseau locations, and a is a vector contain- 
ing the nine coefficients, a.^, of these terms. The quantity 6x is the 
displacement in x at the reseau point, A similar equation yields the y 
displacement coefficients: 

t » *v 

i=0 i=0 

b = A" 1 6y 

These a and b coefficients are computed by the routine entitled SUBREG. 
The BIQUAD routine uses them to compute the displacements for all the 
subreseau points. Corrections for earth curvature and satellite attitude 
are supplied to SUBREG by the INCORP routine. By applying them at this 
level, it is possible to use an adequate second order model without the 
high cost computationally of a later correction to the 4096 2 pixel points. 
Where ground control points are used, the GTRUTH subroutine provides an 
attitude correction vector to INCORP resulting from a combination of at- 
titude and ephemeris predicted data with the ground truth observations. 
When the ground control points are not found or are not available for 
that RBV frame, GTRUTH provides the attitude correction vector based on 
attitude instrumentation and ephemeris data only. 


3.2.4 Bilinear Coefficient 

Using the displacements computed for each subreseau point, sub- 
routine LINCO computes the 8 coefficients (4 for x and 4 for y interpola- 
tion) required for the bilinear distortion model for pseudo reseau 
points. 
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The coefficients are defined by: 


6X, 


'l x, x,y, y, 


V 

ix 2 


1 x 2 x 2 y 2 y 2 


B x 

«x 3 


1 *3 *3*3 y 3 


C x 

6X4 


1 *4 x 4 y 4 y 4 


D x 


The above may be simplified by noting that (x^,y^), {* 2 ^ 2 )' ^ x 3 ,<y 3^’ 
and (x^,y^) are corners of a square with sides of 113. (x^x-i » y 2 = y-| » 
x 2 =x i + H 3, etc.). An analytic inversion of this set of equations for 
the coefficients: 

C x = (l!3) z ^ 6x 4 + 6X 1 " 6x 3 " 6 x 2^ 

B x = TT3 ^ 6x 2 " 6X 1* " C x y l 

D x = TT3 ^ 6x 3 " 6X 1^ " C x X 1 

A x X 1 - B x x l " C x x l y l - D x y l 

For y: C , B , D , A are obtained by analogous expressions with 6x. 

J J J J 

replaced by sy r 


These computations are made in subroutine LINCO. 


3.3 RBV DATA PROCESSING RESULTS 

Due to the shutdown of the RBV subsystem shortly after ERTS-1 
launch, processing priority was shifted to MSS data although considerabl 
progress had been made during the pre-launch study phase to modify TRW 
algorithms to accommodate the RBV data. The software was found to be 
nearly adequate the first time It was utilized and only minor modifi- 
cation was required. 
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With respect to reseau extraction, no measurable position errors 
were found, and no false detections were encountered. Average CPU time 
for reseau extraction was found to be 0.415 seconds, including cases of 
threshold iteration (no more than two such iterations were required for 
any reseau). 

GCP extraction software performance has already been discussed in 
Section 2. Throughput performance for the RBV software using a CDC 6400 
computer is summarized in Table 3-1. 


Table 3-1. Module Throughput 

Element 


CPU Time 

PASS I 



1 . 

Reseaux Extraction 
(81 Reseaux) 

33.6 sec 

2. 

Distortion 

Coefficient Calculation 

10.1 sec 

PASS II 

(Per Band) 


Nearest Neighbor 

8.4 min 

Bilinear 

14.6 min 

TRW Cubic Convolution 

30.8 min 


The image 1002-18134-1 is typical of RBV processing results. Figure 
3.4 illustrates the NASA system corrected image. Figure 3.5 shows a 
detail 960 pixels by 1130 lines reproduced on a line printer from the 
bulk CCTs. Figure 3.6 shows the subregion following precision digital 
processing by means of the methods described earlier in this section. 
Figure 3.7 shows the corresponding detail from the NASA system corrected 
image. Section 4 contains results for various resampling comparisons and 
analyses. 
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$ 


Figure 3.5. Detail from Image 1002-18134-1 Bulk Data 


This figure is a photo reduction of line printer outputs. The 
hash marks running vertically represent position indicators; they 
are not in the processed data. 
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Figure 3.6 Detail from TRW Precision Product for Image 1002-18134-1 
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Figure 3.7. Detail for NASA System Corrected Product for Image 1002-18134-1 
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4 

RESAMPLING 

COMPARISONS 


A significant problem in correcting bulk-image data is the generation 
of image data for pixel locations other than the sampled locations. There 
are two ways to correct image data after determining the necessary geo- 
metric and radiometric corrections: (1) Move the sampled pixels around, 

(e.g., EBR beam position modulation); (2) interpolate the pixel values 
around the desired location ("resampling"). 

TRW has employed the latter process, which in contrast to the former 
produces a digital data stream as the output product. The input to the 
resampling process, which is the second pass process indicated in Figure 1.2, 
consists of the file of distortion coefficients for 64 piecewise bilinear 
distortion regions and the reformatted bulk image data. Conceptually, the 
Pass-II process consists of reconstructing in its entirety the corrected 
(rectified) image and then sampling its intensity values at predesignated 
positions (lines/pixels) within the bulk image. It should be clearly 
understood that at no time is the input image data physically distorted. 

The power of the resampling method becomes clear when it is realized that 
complete flexibility is afforded in the matter of pixel spacings and line 
spacings (independently), thus making possible the custom tailoring of 
the corrected image to any application and/or hardware constraints. 

The resampling method is implemented by means of an interpolation 
process. One-dimensional interpolation of band-limited digital data is of 
the form: 

I(x i> I(u k )f <v u k> 

u k 

wherein f(x) is the interpolation kernel, x. is the argument (pixel loca- 
tion) chosen for the resampled function I, and the set of u k is the set of 
arguments of the available digital data I(u k ), that is, the bulk data 
values. Three kernels, corresponding to: (a) nearest-neighbor interpola- 
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tion; (b) bilinear interpolation; and (c) TRW's Cubic Convolution Pro- 
cess, are shown in Figure 4.1 for one dimension. In the ideal band- 
limited case f(x) is of the form sinx/x, which TRW has approximated by 
the cubic spline function of Figure 4,1-c. 

Nearest-neighbor interpolation involves choosing the value of the 
pixel closest to the desired pixel location as the interpolated value, 
as shown for one dimension in Figure 4.2 . Only one pixel value of the 
bulk image is required for each interpolated value, and only original 
values are used. The algorithm is thus, very fast. Note however that 
errors accumulate in finite (+1/2 pixel) increments and thus will 
result in 1 pixel offsets throughout the corrected image. 

Bilinear interpolation, also shown for one dimension in Figure 4.2, 
utilizes a bilinear combination of the four closest pixel values to pro- 
duce a new, interpolated pixel value. The smoothing effects of bilinear 
interpolation preclude the one pixel offsets characteristic of nearest- 
neighbor interpolation. On the other hand, this smoothing effect causes 
some image degradation in the form of edge smoothing and loss of 
maxima/minima fidelity. Also, the algorithm is inherently slower, inas- 
much as a new pixel value must be computed from four other values. 

In contrast to nearest-neighbor and bilinear interpolation, TRW's 
Cubic Convolution Process requires the values in a grid of 4 pixels x 4 
lines about the point at which the interpolated value is desired. Thus, 
both slope and value continuity properties of the function to be re- 
sampled are preserved (as expected on the basis of its approximation to 
ideal sinx/x resampling). As a consequence, high resolution with minimum 
distortion results, with some penalty in CPU time compared to bilinear 
interpolation. It will be seen, however, that reasonable CPU times have 
been attained by TRW. 

A comparison of performance is given in the error histograms of 
Figure 4.3. A reference image was generated by 900 point sinx/x and 
differenced pixel by pixel with the same image processed by 100 point 
sinx/x, 16 point TRW Cubic Convolution, 4 point bilinear interpolation 
and 1 point nearest-neighbor interpolation. The error histograms shown 
in the figure indicate the clear superiority of TRW Cubic Convolution, 
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if 0 1 j x | i 1 
- | x | if 1 <_ |x 


< 2 


(c) 


Figure 4.1. Three Interpolation Kernels, The three kernels correspond 
to (a) nearestrneighbor interpolation, (b) bilinear 
interpolation, and (c) TRW's Cubic Convolution Process, 
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Figure 4.2. NearestrMeighbor and Bilinear Interpolation 


20634-6003-TU-02 
Page 4-5 


NO. OF 
PIXELS WITH 
GIVEN ERROR 


NO. OF 
PIXELS WITH 
GIVEN ERROR 


lOOOO 


1000 


100 


10 


10000 


1000 


100 


T 

1 

SINC(X+5) 

T 

RMS = 
MAX = 

0.95 

4 

1 









10000 


-1000 


100 


10 


-Jl 


3 


1 1 

TRW CUBIC 
CONVOLUTION - 


□ 

RMS = 0.77 
MAX = 5 


1 




\ 





10 20 0 
-| 1 — 10000 

NEAR-NEIGHBOR 

RMS = 2.62 iooo 
= 26 



10 


20 


\ 

1 

BILINEAR 
RMS = 1.09 


[\ 

MAX = 

9 



\ 




\ 




0 10 20 0 10 20 

ERROR MAGNITUDE ERROR MAGNITUDE 


(INTENSITY RANGE IS 0 ~ 63) 


Figure 4.3. Resampling Errors 
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which generated only one pixel (out of 49,600) with an error greater than 
100 point sinx/x. TRw's Cubic Convolution Process gave the lowest rms er- 
ror of the four considered. 

Examples of imagery produced by the nearest-neighbor, bilinear and 
IRW Cubic Convolution Process are shown in Figure 4.4: the bulk image in 

the upper left, and the processed image using nearest-neighbor interpolation 
(upper right), bilinear interpolation (lower right), and TRW's Cubic Con- 
volution Process (lower left). Note the many one pixel image discontinu 
ities characteristic of nearest-neighbor interpolation, particularly evi- 
dent for the road intersection in the upper left corner of the processed 
image. Bilinear interpolation, on the other hand, eliminates these dis- 
continuities, at the expense of image resolution. Finally, the image pro- 
cessed by the Cubic convolution Process shows none of the nearest-neighbor 
image discontinuities, and no loss of resolution. Figures 4.5 and 4.6 more 
clearly show the differences between nearest-neighbor and the TRW Cubic 
Convolution Process, respectively, reproduced by a filmwriter. Figure 4.7 
shows the bulk image equivalent for the same area, derived from bulk CCT 
data. 

Spatial frequency analyses based upon the power density spectra for 
processed ERTS data have also been performed. Figure 4.8 was chosen as a 
typical image detail (scene 1062-15190-5), 100 lines x l 28 pixels in the 
corrected image. This detail shows the intersection of 183 (Jones Falls 
Expressway) and 1695 (Baltimore Beltway), and is reproduced from line 
printer output.* Using TRW's Fast Fourier Transform Process, line power 
density spectra were obtained for the detail and then averaged (spectral 
component by spectral component) over all 100 lines. Tnis was done in turn 
for nearest-neighbor, oil inear and TRW's CuDic Convolution Process. To 
eliminate leakage due to a large d- component, the mean was removed, pad- 
ding zeros were added to the data and a Hamming Window (cosine taper at tne 
image edges) was employed. 

The results are plotted as a function of |w| (rad/meter) in Fig- 
ure 4.9. The low frequency contents of three resampled images are very 

^ " 1 

A dot matrix 3 dots on a side permits use of 10 gray levels for each 
pixel . 
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Figure 4.4. Corrected Image Detail for Three Interpolation Algorithms 
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Figure 4.5. Nearest-Neighbor Processed Image Detail 
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Figure 4.6. TRW Cubic Convolution Processed Image Detail 
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Figure 4-7. BULK IMAGE SUBAREA FROM 1062-15190-5. 
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Figure 4,8. Detail of Precision Processed 1062^1 51 90t.5 

This 1 iner.printer output shows the Intersection of 183 and 1695, 
resulting from TRW Cubic Convolution' 'Processing. 



NEAREST NEIGHBOR 


Figure 4. 


. POWER SPECTRAL DENSITY: AVERAGED FOR 100 LINES, ERTS WARP 

The folding frequency Is denoted by and the 
frequency Interval fa => radian/meter 
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nearly the same. (However, that corresponding to the TRW Cubic Convolution 
Process contains more high frequency components than bilinear resampled 
images. This is an expected result due to the smoothing nature of bi- 
linear interpolation. The nearest-neighbor processed image is not greatly 
different from the TRW Cubic Convolution (bear in mind, the small relative 
power at the highest spatial frequencies). 

Another way of comparing the three resampling techniques is to 
generate three images shifted one-half pixel for the scene of Figure 4.10 
(bulk data) by the use of nearest-neighbor, bilinear and TRW Cubic 
Convolution Process, respectively. The power density spectra of the three 
shifted images were computed by again averaging over 100 lines. Results, 
plotted in Figure 4.11 , show that the power density spectrum of nearest- 
neighbor resampled images is indistinguishable from the power spectrum of 
the original bulk image since a linear shift of half pixel by means of 
nearest-neighbor resampling should preserve the image samples and thus the 
statistics. The bilinear resampled image suffers the loss of high frequency 
components because of the smoothing effect of bilinear reconstruction. The 
power density spectrum corresponding to TRW's Cubic Convolution Process has 
essentially the same spatial frequency content as the nearest-neighbor 
resampled image; however, the high frequency components are not degraded 
as is the case for bilinear interpolation. 
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igure 4.10. Second Detail of Precision Processed 1062-15190-5 



Relative Power of Image 
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